In [1]:
# Setting up a model and a mesh for the MT forward problem
In [5]:
import SimPEG as simpeg, sys
import numpy as np
from SimPEG import NSEM
sys.path.append('/Volumes/MacintoshHD/Users/thibautastic/Dropbox/PhD_UBC/telluricpy')
import telluricpy
import matplotlib.pyplot as plt
%matplotlib inline
import copy
In [6]:
# Define the area of interest
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480
In [7]:
# Build the mesh
# Design the tensors
hSize,vSize = 25., 10
nrCcore = [15, 8, 6, 5, 4, 2, 2, 2, 2]
hPad = simpeg.Utils.meshTensor([(hSize,10,1.5)])
hx = np.concatenate((hPad[::-1],np.ones(((be-bw)/hSize,))*hSize,hPad))
hy = np.concatenate((hPad[::-1],np.ones(((bn-bs)/hSize,))*hSize,hPad))
airPad = simpeg.Utils.meshTensor([(vSize,13,1.5)])
vCore = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcore,(simpeg.Utils.meshTensor([(vSize,1),(vSize,8,1.3)])))])[::-1]
botPad = simpeg.Utils.meshTensor([(vCore[0],8,-1.5)])
hz = np.concatenate((botPad,vCore,airPad))
# Calculate the x0 point
x0 = np.array([bw-np.sum(hPad),bs-np.sum(hPad),bt-np.sum(vCore)-np.sum(botPad)])
# Make the mesh
meshFor = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
In [8]:
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
meshFor=copy.deepcopy(meshInv)
In [9]:
print np.sum(vCore)
print meshFor.nC
print meshFor
In [23]:
# Save the mesh
#meshFor.writeVTK('nsmesh_FineHKPK1.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
In [24]:
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)
In [25]:
#telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)
In [26]:
# Get active indieces
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')
In [27]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
'PK1':5e-2,
'HK1':1e-3,
'VK':5e-3}
In [28]:
# Loop through
extP = '../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')
In [29]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','PK1','HK1','VK']:
sigma[geoStructIndDict[key]] = geoStructFileDict[key]
In [37]:
fig = plt.figure(figsize=(10,8))
ax = plt.subplot(111)
model = sigma.reshape(meshFor.vnC,order='F')
a = ax.pcolormesh(meshFor.gridCC[:,0].reshape(meshFor.vnC,order='F')[:,20,:],meshFor.gridCC[:,2].reshape(meshFor.vnC,order='F')[:,20,:],np.log10(model[:,20,:]),edgecolor='k')
ax.set_xlim([bw,be])
ax.set_ylim([0,bt])
ax.grid(which="major")
cb = plt.colorbar(a)
ax.set_aspect("equal")
cb.set_label("Log conductivity (S/m)",fontsize=20)
ax.set_xlabel("Easting (m)",fontsize=20)
ax.set_ylabel("Elevation (m)",fontsize=20)
plt.tight_layout()
$\nabla ^2 \mathbf{H} + k^2 \mathbf{H} = 0$
In [38]:
# Save the model
meshFor.writeVTK('nsmesh_FineHKPK1Fin.vtr',{'S/m':sigma})
In [39]:
# Set up the forward modeling
freq = np.logspace(5,2,16)
np.save('MTfrequencies',freq)
In [50]:
telluricpy.vtkTools.extraction.geometryFilt??
In [40]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
import vtk
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw-5.,be+5,bs-5.,bn+5.,bb-5.,bt+5)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=False)
telluricpy.vtkTools.io.writeVTPFile('topoSurf.vtp',actModVTP)
In [52]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+10.,be,20),np.arange(bs+10.,bn,20))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,actModVTP) #modTopoVTP)
np.save('MTlocations',locArr)
In [53]:
telluricpy.vtkTools.io.writeVTPFile('MTloc.vtp',telluricpy.dataFiles.XYZtools.makeCylinderPtsVTP(locArr,5,10,10))
In [60]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py
In [71]:
%matplotlib qt
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf
In [72]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_HKPK1.npy')
mtData
Out[72]:
In [110]:
iPf.MTinteractiveMap([mtData])
Out[110]:
In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions,
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.
Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py
In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')
In [109]:
np.unique(drecAll['freq'])[10::]
Out[109]:
In [93]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
In [94]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC)})
In [101]:
nsInvvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
activeModInv = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsInvvtr,topoSurf,extBoundaryCells=True)
In [102]:
sigma = np.ones(meshInv.nC)*1e-8
indAct = telluricpy.vtkTools.dataset.getDataArray(activeModInv,'id')
sigma[indAct] = 1e-4
In [103]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC),'S/m':sigma})
In [108]:
from pymatsolver import MumpsSolver
In [106]:
pymatsolver.AvailableSolvers
Out[106]:
In [116]:
NSEM.Utils.skindepth(1000,100000)
Out[116]:
In [117]:
np.unique(mtData['freq'])[5::]
Out[117]:
In [ ]: